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An iterative method is derived for image reconstruction. Among other attributes, this method 
allows constraints unrelated to the radiation measurements to be incorporated into the reconstructed 
image. A comparison is made with the widely used Maximum-Likelihood Expectation-Maximization 
(MLEM) algorithm. 
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Imaging by radiation emission or transmission effec- 
tively produces a set of linear equations to be solved. For 
example, in the case of coded aperture imaging, the so- 
lution is a "reconstructed" set of radiation sources, while 
in the case of x-ray interrogation, the solution is a set 
of attenuation coefficients for the voxels comprising the 
volume through which the x-ray beam passes. 

The linear equations have the form 



d i = Y d M i 
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where the set {di} corresponds to the radiation intensity 
distribution recorded at a detector (a detector pixel is 
labeled by the index i), the set {^j} is the solution, and 
the matrix element My connects the known di to the un- 
known . Typically the matrix M is non-square so that 
{fJ,j} cannot be obtained by standard matrix methods. 
(And note that, when the set of equations is large, it can 
be difficult to ascertain a priori whether the equation set 
is over- or under-determined.) 

In any case the matrix equation d = Mfi may be 
solved by the iterative method that is derived as follows. 
Clearly this method will feature a relation between 

and fij , where n is the iteration number. Consider 

the two equations for [a^ and fj/^ 1 ] 



and rewrite the latter as 
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Then the relationship between fj^ and (Aj ^ is ob- 
tained by setting X^df^ = Ei<^ : 
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Note that this last equation can be written 



(5) 
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where the last factor is essentially a weighted average of 
all di/dj Tl ~ 1 \ Thus the set i/^™ j a PP r °aches a solution 

{/ij} by requiring Y]jd^ = Y)jdj at each iteration; in 
effect, by requiring all — > di. 

The iteration procedure alternates between use of Eq. 
© and Eq. ([5]) until all g^™' are as close to di as desired. 
For the first (n = 1) iteration, an initial set {mj ' j is 

chosen, which produces the set {di "*} according to Eq. 
©. These values are used in Eq. ([5]), so producing 
the set j/x^|. And so on... That a final set |/Zj- n H 
is a solution {fj,j} to the matrix equation d = M/j, is 
verified by checking that all d\ = di to within a desired 
tolerance. 

Some cautions and opportunities follow from this sim- 
ple derivation of Eq. ©. A caution is that, in the event 
the equation set is under-determined, different initial sets 
j/u^ j will lead to different final sets {/Xj} that satisfy 
the matrix equation. The corresponding opportunity is 
that this problem may be mitigated to some extent by 
the addition, to the original set of equations, of linear 
equations that further constrain the (perhaps derived 
from, for example, independent knowledge of some of the 
contents of a container under interrogation). In general 
the di appearing in a constraint equation will have noth- 
ing to do with radiation intensity. 

The form of any added constraints, and the initial 
choice {mj " 1 }' mus t allow all ~ * Mj an d <4™' — > di 
monotonically. In particular, care should be taken when 
a constraint has one or more coefficients My < 0, as that 
affects the denominator Yli-^-ij i n Eq. © (a straightfor- 
ward fix may be to reduce the magnitudes of all My 
coefficients and di in that constraint equation by a mul- 
tiplicative factor). In any event, the acceptability of a 
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constraint equation is easily ascertained by monitoring 

(n) 

the behavior d\ —> di for that constraint. 

Note that all solutions {/ij} to a set of equations 
that includes additional constraints with di > and all 
Mij > are accessible from sets {mj ' j °f initial values, 

and further that any set | /i^ "* | will produce a solution 

{fJ-j}- This suggests that, for this implementation of con- 
straints, a superposition of many solutions may give a 
good "probabilistic" reconstruction. To achieve this, con- 
sider that the innumerable solutions to the set of equa- 
tions may be regarded as points in a J-dimensional space 
(J is the number of elements in a solution {/%}). These 
points must more-or-less cluster, producing a cluster cen- 
troid that is itself a solution. While the centroid solution 
{/^ J has no intrinsic special status (as all cluster points 

represent equally likely reconstructions) , it may be taken 
to represent the particular set of equations. The cluster 
size, which indicates the degree to which solutions are 
similar to one other, should decrease as constraints are 
added. A logical measure of the cluster size is 

^cluster = ((x c - Xfc) • (x c - Xfc))^ 2 (7) 

where x c is the centroid vector and x^ is the vector 
corresponding to the fcth solution. Thus the quantity 
fciustcr/v^; which represents the standard deviation of 
the innumerable values of an arbitrary element fj,j, is a 
useful measure of the variation among solutions {fJ,j}. In 
general it is desirable that the variation among solutions 
be much less than the variation within the centroid solu- 
tion, which is 




where ^( c ) = (Mj / .■ I R that case (^cluster-/ -C ojf ) 
the centroid solution is little changed by additional con- 
straints, so suggesting that the centroid solution j 1 



may be regarded as the sought-after reconstruction. 

Another caution follows from the fact that the denom- 
inator J^i^ij in Eq. ([5]) is the sum of all elements in col- 
umn j of matrix M. This iterative method can of course 
be used without explicitly converting a set of linear equa- 
tions into a matrix equation (or several sets of equations 
into a single matrix equation), but in that event very 
careful attention must be paid to get the factors Yli-^ij 
right. 

It may be noticed that Eq. ([5]) is similar to the so- 
called Maximum-Likelihood Expectation-Maximization 
(MLEM) algorithm (see refs. [l[ and for derivations of 
the latter, and see numerous papers in the recent imag- 
ing literature for applications of it). The MLEM, which 
is derived from physical considerations having to do with 
radiation emission and detection, purports to find the 
set {nj} that maximizes the probability P ({di} | {/%}), 
which is the probability of realizing the observed set {di} 
given a set {nj}- This is in contrast to the derivation 
above, which leads to Eq. ([5]) as simply a method to find 
a solution to a matrix equation (a set of linear equations). 
The derivation presented here makes clear how the itera- 
tive procedure should be implemented for an application, 
and allows constraints to be added to the original set 
of equations (those produced by the imaging exercise) 
thereby enabling a more-accurate reconstruction when 
the original equation set is under-determined. 
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